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The question of reliability arises for any dynamical system driven by an input signal: if the same 
signal is presented many times with different initial conditions, will the system entrain to the signal 
in a repeatable way? Reliability is of particular interest in large, randomly coupled networks of 
excitatory and inhibitory units. Such networks are ubiquitous in neuroscience, but are known to 
autonomously produce strongly chaotic dynamics - an obvious threat to reliability. Here, we show 
that such chaos also occurs in the presence of weak and strong stimuli. However, even in the chaotic 
regime, intermittent periods of highly reliable spiking often coexist with unreliable activity. We 
argue that the sustained coexistence of chaos and reliable spike events is due to the interaction of 
global state space expansion and dynamics local to individual cells, and interpret our findings within 
the framework of random dynamical systems theory. 



Reliability is a general property of dynamical systems 
driven by input signals. This refers to the reproducibility 
of a system's output on many trials, where the same driv- 
ing signal is presented but initial system states vary. The 
concept is fundamental in neuroscience: for a network of 
interacting cells receiving sensory and internal stimuli, 
the degree to which it is reliable impacts its ability to 
encode information via precise temporal spike patterns. 
Related phenomena arise in a variety of physical and en- 
gineered systems, including coupled lasers pQ and cou- 
pled chaotic systems [2]. Understanding the conditions 
and dynamical mechanisms that govern network reliabil- 
ity stands as a broad challenge in the study of networks 
of dynamical systems. 

The question of reliability is of particular relevance for 
an ubiquitous and important class of neural networks, 
those with a balance of excitatory and inhibitory connec- 
tions [3J . Such balanced networks produce dynamics that 
match the irregular firing observed experimentally on the 
"microscale" of single cells, and on the macroscale have 
rapid and linear mean-field dynamics that could be ben- 
eficial for neural computation 0HH]- However, such bal- 
anced networks are known to be produce strongly chaotic 
activity when they fire autonomously or with constant in- 
puts [6l0[9]. On the surface, this appears incompatible 
with reliable spiking, as small differences in initial condi- 
tions between trials may lead to very different responses. 
However, that the answer might be more subtle is sug- 
gested by a variety of results on the impact of temporally 
fluctuating inputs on chaotic dynamics [T0l - fT6] . 

In this Letter, we present a detailed numerical study 
and steps toward a qualitative theory of reliability in 
fluctuation-driven balanced networks. We explore the 
relationship between the Lyapunov spectrum, which 
quantifies the average stability of trajectories on long 
timescales, and the cross-trial repeatability of spike times 
over shorter timescales. Even in the presence of positive 
Lyapunov exponents, we find that spike times can dis- 
play sharp temporal precision from trial to trial. These 
reliable spike events are interspersed with more diffuse 
spiking, in which spike times are impacted by prior net- 



work states differing from trial to trial. Thus, networks 
with positive Lyapunov exponents produce two types of 
spike events, representing either inputs alone or a com- 
bination of inputs and past network states. 
Model and single-trial dynamics: We study a tempo- 
rally driven network of N — 1000 spiking neurons. Each 
neuron is modeled by a phase variable 9i € S 1 — R/Z 
whose dynamics follow the "9- neuron" model [T7], cap- 
turing the saddle-node on an invariant circle (SNIC) bi- 
furcation responsible for spike generation in Type I neu- 
rons. These dynamics are equivalent to the quadratic 
integrate-and-fire (QIF) model after a change of coordi- 
nates jTHj. The underlying "normal form" dynamics [TS] 
are found in many brain areas and are known to pro- 
duce reliable responses to stimuli in isolated cells [T5ll20| . 
cf. [31] [32]. Thus, any unreliability or chaos that we find 
is purely a consequence of network interactions. 

Coupling from neuron j to neuron % is determined by 
the weight matrix A = {dy}. A is chosen randomly using 
an Erdos-Renyi scheme such that 20% of the cells j are 
inhibitory (<Zjj < Vi) and 80% are excitatory (ay > 
Vi); we do not allow self-connections, setting an = 
0. Each neuron has mean in-degree K = 20 from each 
population (excitatory and inhibitory) and the synaptic 
weights are 0(l/y/~K) in accordance with the classical 
balanced-state network architecture [6]. We note that 
our results appear to be qualitatively robust to changes 
in N and K, but a detailed study of scaling limits is 
beyond the scope of this paper. 

A neuron j is said to fire a spike when 6j(t) crosses 
9j = 1; when this occurs, 9i is impacted via the coupling 
term atjg(9j) where g{9) is a smooth "bump" function 
with small support ([—1/20,1/20]) around 9 = sat- 
isfying J g{9)d9 = 1, meant to model the rapid rise 
and fall of a pulsatile synaptic variable. In addition 
to coupling interactions, each cell receives a stimulus 
h (t) = i] + eQ (t) where 77 represents a constant current 
and Q(t) are independent white noise processes, scaled 
by the amplitude parameter e. Here, C,i{t) represents a 
"frozen" (quenched), aperiodic signal (as in [T3I |2"TI |2"2"] ) . 

The z th neuron in the network is therefore described 
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by the following SDE: 

Mi = [F(6i) + Z{0i) \V + Y1 a H9( e j) + 

V i J (1) 

-z{6 i )z'{e i )]dt + ez{e i )dw i ,t 

where the intrinsic dynamics F(6i) — 1 + cos(27rf9j) and 
the stimulus response curve Z(0i) = 1 — cos(27rf3i) come 
directly from coordinate changes from the original QIF 
equations (see Supplemental Materials (SM) [23] and 
|17j). Here, Wi t is the independent Wiener process gen- 
erating d(t); the e 2 term is the Ito correction from the 
coordinate change [24]. Finally, r\ sets the intrinsic ex- 
citability of individual cells. For i] < 0, there is a sta- 
ble and an unstable fixed point, together representing 
resting and threshold potentials. Thus (contrasting [8 
where cells are intrinsically oscillatory), cells are in the 
"excitable regime," displaying fluctuation-driven firing, 
as for many cortical neurons |25j . 

Fig. [T] illustrates that the general properties of the net- 
work dynamics, including a wide distribution of firing 
rates from cell to cell and highly irregular firing in in- 
dividual cells, are consistent with many balanced-state 
networks from the literature as well as general observa- 
tions from cortex [H[S]. Furthermore, the network's mean 
firing rate scales monotonically with rj and e (not shown) 
as in [H] . 

Asymptotic reliability: We say a network is asymp- 
totically reliable if for any fixed input signal (a realization 
of Q(t) = {Ci(i)}i=i...iv)j solutions starting with distinct 
initial conditions (ICs) converge to a single trajectory, 



A Typical firing rate histogram B Typical single-cell ISI distr. 




FIG. 1. (A) Typical firing rate distributions for excitatory 
and inhibitory populations. (B) Typical inter-spike-interval 
(ISI) distribution of a single cell. The coefficient of Variation 
(CV) is close to 1. (C) Invariant measure for an excitable 
cell (n < 0); inset: typical trajectory trace of an excitable 
cell where solid and dotted lines mark the stable and unsta- 
ble fixed points. (D) Network raster plots for 250 randomly 
chosen cells. For all panels, n — — .5, e = 0.25. 



therefore producing spikes at the same precise time on 
every trial. 

In this context, it is useful to treat ([I]) as a random 
dynamical system (RDS) on T N . That is, we view the 
system as a nonautonomous ODE driven by a frozen re- 
alization of £(t), and consider the action of the generated 
family of flow maps on phase space. For such a RDS, Lya- 
punov exponents Ai > A2 > ... > \n are well defined and 
are constant for almost every choice of f(i) and almost 
every IC. Assuming a number of nondegeneracy condi- 
tions, it can be shown [5S1 [57] that when Ai < 0, almost 
every trajectory will converge to a random sink: a in- 
dependent, asymptotically stable trajectory on T N . We 
note that although deterministic systems tend to have 
multiple basins with Ai < 0, RDS's often have a single 
random sink because of ergodicity (see SM [13]). On the 
other hand if Ai > 0, solutions are attracted to random 
strange attractors [2S] (called "snapshot attractors" in 
[29]). These are time-dependent versions of Sinai- Ruelle- 
Bowen measures for dissipative chaotic systems |30] , and 
often are not localized in state space. We therefore use 
the sign of Ai as an indicator of the asymptotic reliability 
of the network ([I]). That is, if Ai < 0, we expect almost 
every solution to converge to a unique trajectory whereas 
if Ai > 0, this is a strong indication that dynamics are 
chaotic and thus not asymptotically reliable. 

In Fig. [2] (A) we present the dependence of the first 100 
Lyapunov exponents for a network with fixed mean in- 
put 77 = —0.5 but varying fluctuation amplitude £ (see 
SM [23] for numerical details). Fig. [2] (C) shows the 
dependence of Ai on e. The networks produce a neg- 
ative Ai for sufficiently small values of the mean input 
e. However, as e increases Ai becomes positive, indi- 
cating chaotic network dynamics and thus, asymptotic 
unreliability. To better understand the nature of this 
chaotic regime, we compute the Kaplan-Yorke dimension 
(Dky — k+Si=i -Wl^fe+il where k is the largest integer 
such that 53i=i Ai > 0), which characterizes the strange 
attractor on which the dynamics evolve. Fig. [2] shows 
Dky/N as a function of e, representing the fraction of 
phase space on which the dynamics are concentrated. In- 
terestingly, even as Ai increases with e, D KY /N eventu- 
ally decays as e gets bigger. 

We now explore the effects on network outputs. 
Spike time reliability: We use a metric i? S pike which 
quantifies the similarity of a family of spike trains by 
measuring the fraction of spikes that are repeated in 
each train. If all spike trains arc identical, i? S pikc — f 
while if all spikes are far enough apart in time, i? S pike ap- 
proaches zero. To compute this value, we follow [3T] and 
define cross-trial spike events by time windows surround- 
ing peaks in the spike times histogram, obtained from 
combining filtered spike trains (see raster plots in Fig. [2] 
(B) and SM [23]). Events that contain a spike from each 
train are labeled reliable, as are the spikes in them. 

Using this metric, we compare the spike output of in- 
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FIG. 2. (A) First 100 Lyapunov exponents of network with 
parameters as in Fig. [TJ as a function of e. (B) Raster plots 
show example spike times of an arbitrarily chosen cell in the 
network on 30 distinct trials, initialized with random ICs. 
Circle and star markers indicate e values of 0.18 and 0.45, 
respectively, shown in panel (C). (C) Multi-scaled plot of Ai 
(right scale), Dky/N and 1 — i^p^e (left scale) vs e. Error 
bars on 1 — i? ap ike curve indicate standard deviation of value 
across all cells in the network. (D) IRI log histogram (see 
text) for a sample parameter set (star marker). (E) Variance 
vs mean for cross-trial single cell spike counts in a time win- 
dow of width T=2 (similar results for different time windows) 
for a sample parameter set (star marker). Dotted line shows 
identity (unit Fano factor). For all panels, r) = —0.5. 



dividual cells in the network on 30 trials initialized with 
randomly chosen, network-wide ICs but receiving the 
same input set Q(t). Panel (C) of Fig. [2] shows 1 — i? S pike 
versus e, averaged over every cell in the network. We 
compute this curve using 2500 time unit runs, discarding 
the initial 10% to avoid transient effects. As expected, 
parameter regions where Ai < correspond to i? sp ikc = 1 
while if Ai > (the chaotic regime), i? sp ikc < 1- 

However, even if the network is chaotic, i? S pikc remains 
positive, indicating that many spikes are still reliably 
repeated across trials. Panel (B) of Fig. [2] illustrates 
this via raster plots of cross-trial spike times of a ran- 
domly chosen cell: even if Ai > 0, the spiking shows 
repeated temporal structure from trial to trial. The tem- 
poral statistics of reliable spike events follow an expo- 
nential distribution, as pictured in Fig. [2jD) showing 
the log histogram of inter-reliable event intervals (IRI). 
Furthermore, this phenomenon induces cross-trial spike 
count statistics that further suggest temporal precision: 
Fig. [2|E) compares the spike count variance and mean, 
computed across trials in time windows of length T = 2. 
When Xi < 0, the spike count variance is zero as ex- 
pected (not shown). When Ai > 0, the variance remains 
generally lower than the mean. This indicates a Fano 



factor less than one, thus suggesting spike time precision 
that substantially exceeds that of the benchmark (inho- 
mogeneous) Poisson process [32] . 

Importantly, we note that the reliable spike events we 
observe do not simply follow from choosing the external 
input amplitude e to completely dominate all network in- 
teractions. Rather, for all parameters explored, network 
interactions generate larger fluctuations in cell responses 
than the external input eQ(t) (see SM 23J). Thus the co- 
existence of chaos and reliable spike events is a dynamical 
phenomenon, whose origins we now investigate in more 
detail. 

Coexistence of reliable spike events and chaos 

for Ai > 0: The persistence of reliable spike events in 
chaotic regimes may be surprising, as one generally ex- 
pects trajectories to diverge at an exponential rate Ai. 
However, this divergence is confined to lower-dimensional 
subspaces characterized by the random strange attractor, 
and is heterogeneous in both time and "space" . We sug- 
gest this explains in part the trend in Fig. [2] (B), where 
1 — i? sp ike more closely resembles Dky/N than Ai. 

To better understand the effect of localized space ex- 
pansion on shorter timescales, we select an example pa- 
rameter set for asymptotically unreliable dynamics - 
listed in the caption of Fig. [3]- and show that spikes in 
unreliable events are associated with a greater degree of 
expansion in the direction of the spiking cells. Consider 
v(t), the solution of the variational equation v = J(t)v 
where J(t) is the Jacobian of the flow evaluated along the 
trajectory 9{t). If we set i>(0) to be randomly chosen but 
with unit length, then v{t) quickly aligns to the directions 
of maximal expansion in the tangent space of the flow 
about 9(t); moreover, Ai = limbec \ log(||u(t)||). We can 
equivalently write a discretized version of this expression 
for small At: Ai = lim/r_ >tx) (e(i)}T where (-)t denotes the 
time average up to time T and e(t) = ^ log 
is analogous to a finite time Lyapunov exponent. For our 
network, e(i) fluctuates rapidly and depends on many 
factors such as number of spikes fired, the pattern of the 
inputs, and the phase coordinate of each cell over the 
time At. Its coefficient of variation is typically O(10) 
for At = 0.005 which shows that stability is very het- 
erogeneous in time. This further suggests that Ai only 
captures the asymptotic mean of a very wide distribu- 
tion, giving us little insight into the behavior of the flow 
on finite timescales. 

To better understand the impact of the flow on a single 
cell's subspace, we define two quantities 
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The support score Si(t) measures the normalized contri- 
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FIG. 3. For (A-C), t = marks the spike time and rel/unrel 
indicate the identity of the spike used in the average. (A) 
Spike triggered Lyapunov vector support S. (B) Spike trig- 
gered local expansion measure E. (C) Spike triggered average 
phases with unstable fixed point U marking threshold. (D) 
Left: two snapshots of Si(t) for different times. Right: his- 
togram of Si (t) for a randomly chosen cell over a i = 200 long 
trajectory. (E) and (F) Time evolution of distance between 
two distinct trajectories f? 1 (t), 6 2 (t) (E) Green dashed (bot- 
tom): \\6j(t) — 8i(t)\\ s i in randomly chosen 9i subspace. Black 
solid (top): m^i{\\el{t) - ef{t)\\ s i}. (F) p 1 (t) - 2 (t)|| T « . 
Network parameters: c = 0, 7} = —0.5, e = 0.25 with A « 1.5. 



bution of a single cell's subspace to the support of the 
maximal Lyapunov vector v(t). The local expansion co- 
efficient ei (t) is a local equivalent of e(t) . We compare the 
time course of these two quantities in the moments pre- 
ceding spikes, by defining the corresponding spike trig- 
gered averages E(t) and S(t), i.e., the means of e^(t) 
and Si (t) in a short time interval preceding each spike in 
the network. We separate reliable spikes from unreliable 
ones to obtain two versions of both measurements. Fig. [3] 
shows the resulting averages. Moments before a cell fires 
an unreliable spike, S(t) is considerably larger than in 
the reliable spike case, thus indicating that global expan- 
sion is further localized in unreliable spike events. Addi- 
tionally, Fig. [3] (B) shows that J5 un rei W's peak is much 
broader than E le i(t) 's with f_ 2 E unrc i(t) — E ro \(t)dt ~ 2.5 
which indicates that prior to an unreliable spike, trajec- 
tories are subject to an accumulated infinitesimal expan- 
sion rate higher than in the reliable spike case. 

The source of "local" expansion e,(i) is closely linked 
the the phase trajectory 6i(t) prior to a spike. If Gift) < 
|, the derivative of the single cell flow defined in (Tip -in 
absence of fluctuating inputs (from network or external 
sources)- is negative, and becomes positive for 6i{t) > |. 
When an uncoupled cell is driven by Q, we know that 
on average, it spends more time in its contractive region 
(6i < i) and is reliable as a result [T31 [5D] . While inputs 
may directly contribute to J(t), their effect is generally 
so brief that their chief contribution to ei(t) is to steer 



0i(t) in expansive regions of its own subspace. Fig.[3](C) 
confirms that the average phase of a cell preceding an 
unreliable spike spends more time in its expansive region. 
Such a phenomenon has previously been reported in the 
form of a threshold crossing velocity argument [33]. We 
refer the reader to the SM [53J for a detailed treatment 
of input conditions leading to this phenomenon. 

The key feature of this system, likely due to sparse and 
rapid coupling, is that Ai represents a sustained balance 
between inputs leading to contraction/expansion in local 
neural subspaces. A bias toward higher frequencies of 
"expansive inputs" yields Ai > and implies on average 
more growth than decay in the maximal Lyapunov direc- 
tion v(t). Because of the local nature of expansive events, 
this direction is only supported on a few neural directions 
at a time, as reported in [H] and illustrated in Fig. [3] (D), 
where we plot snapshots of u(£)'s components {si(i)}j. 
Furthermore, any given cell seldom contributes to v(t) 
(see Si histogram in Fig. [3](D)) which indicates that the 
identity of subspaces involved in v(t) changes constantly 
(as also in [5]). This leads to trajectories that are un- 
stable on long timescales (Ai > 0), yet alternate between 
periods of stability and instability on finite timescales. 
Fig. [3] (E) shows the time trace of \\6}(t) - 9f(t)\\ s i: 
the projection distance between two randomly initialized 
trajectories and 2 (t) in a single cell's direction as 

well as mDXi{\\6j(t) - 6f(t)\\ s i}. While the maximal S 1 
distance is almost always close to its maximum 0.5, the 
two trajectories regularly collapse arbitrarily close in any 
given S^-subspace. In essence, neurons take turns con- 
tributing to global separation of trajectories but can be- 
come locally stable otherwise. This leads to a global sep- 
aration ||6 |1 (t) — f5 2 (i)|| T ]v that is relatively stable in time 
(Fig.j3j(F)) yet supports the local, sustained coexistence 
of reliable and unreliable spike events. 

Conclusion: In this Letter, we have explored the re- 
liability of fluctuation-driven networks in the excitable 
regime -where single cell dynamics contain stable fixed 
points. Unreliable spikes are a hallmark of sensitivity 
to initial conditions and may therefore carry information 
about previous states of the system (or, equivalently, pre- 
vious inputs). On the other hand, reliable spikes carry 
repeatable information about the external stimulus I(t) 
alone. We showed that both unreliable and reliable spike 
events coexist in chaotic regimes of the system explored. 
The resulting implications for the neural encoding of sig- 
nals are an intriguing avenue for future work. 
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